Efficient approach for simulating distorted materials 
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The operation principles of nanoscale devices are based upon both electronic and mechanical 
properties of materials. Because these properties can be coupled, they need to be investigated si- 
multaneously. At this moment, however, the electronic structure calculations with custom-made 
long-range mechanical distortions are impossible, or expensive at best. Here we present a unified 
formalism to solve exactly the electronic structures of nanomaterials with versatile distortions. We 
illustrate the formalism by investigating twisted armchair graphene nanoribbons with the least pos- 
sible number of atoms. Apart from enabling versatile material distortions, the formalism is capable 
of reducing computational costs orders of magnitude in various areas of science and engineering. 
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Bloch's theorem has been the propulsive force of com- 
putational materials research for more than 80 years [1, 
today as important as ever. While the theorem still asso- 
ciates with Bravais lattices and translational symmetry, 
nanoscience has brought us low-dimensional structures, 
tubes, tori, wires, and membranes, which get twisted, 
bent, wrapped, and rippled in experiments. Transla- 
tional symmetry hides deceiving simulation constraints, 
since materials cannot distort the way they would prefer, 
and restricts realistic modeling of nanoelectromechanical 
components. 

Distortions are relevant in a number of topical mate- 
rial systems: polymers, double helices like DNA, lipid 
bilayers, nanoscrolls, nanocoils, nanowires, and, espe- 
cially, carbon nanostructures including fullerenes, carbon 
nanotubes (CNTs), graphene, and graphene nanoribbons 
(GNRs), to mention a few.^] For example, materials with 
high aspect ratio like CNTs and GNRs get bent [3,-^ and 
thin sheets like graphene get rippled [3 [H [7] , unless care- 
fully placed on a support. Classical modeling of distor- 
tions is a mature subject [8 10 , but while classical inter- 
action potentials and finite element methods give materi- 
als' mechanical properties, they are useless for electronic 
properties. In nanoscience quantum- mechanical model- 
ing is preferred. 

How can we include quantum mechanics into these 
distortion simulations? For decades chemists have used 
group theory and molecular symmetries to reduce compu- 
tational costs. In computational materials physics, sym- 
metries beyond translation have been used mainly for 
chiral carbon nanotubes, in work pioneered by White, 
Robertson and Mintmire[TT], followed by Popov [12] and 
Dumitricap!3], with co-workers. Nanotubes are natural 
because chiral symmetry itself suggests "symmetry adap- 
tion" ; it is, however, less evident to break the symmetry 
and investigate the elastic properties in a broader sense. 

In this Letter we shall present a compact, exact, and 
flexible formalism to solve the electronic structure of 
nanomaterials with custom-made distortions. By ex- 
panding the concepts of periodicity and simulation cells, 
the formalism can also reduce computational costs, even 
in classical materials modeling. 
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The formalism is obtained by revising Bloch's theorem, 
and the derivation is straightforward. Consider electrons 
in a potential V{r) that remains invariant in symmetry 
operations S'^^ 

Z)(5^)V(r) = V(5-^r) = V{r). (1) 

The operation 5^^, with inverse 5"^^, is a succession 
of rii times operation Si for all i = 1,2,..., that is 
S'^ = 81^82'^ ' • ' with n = (ni, n2, • • • )• can be any 
symmetry operation, such as translation, rotation, re- 
flection, inversion, joined translation+rotation, or joined 
translation+reflection, to mention six, and they should 
form an abelian group. 

Let us give a couple of familiar examples. With bulk 
tSi's are three translations; for a benzene ring {CqHg) Si 
could be a two-, three-, or six- fold rotation around the 
symmetry axis; for an achiral carbon nanotube Si could 
be a translation along the symmetry axis and ^2 could 
be, say, an M-fold rotation around the symmetry axis; 
for polyethene ([-CH2CH2-]n) Si could be a translation 
across one CH2 unit followed by a reflection (normally Si 
would be a translation across the whole CH2CII2 unit). 

Now, returning to the derivation, since the transforma- 
tions are isometric, the kinetic energy term in the Hamil- 
tonian 

H = -^V^ + Vir) (2) 

remains invariant, and H commutes with D{S'^)^ the two 
operators consequently sharing the same eigenstates. We 
denote these eigenstates ijjanir)^ with k = (a^i, ^^2, • • •)• 
Hence we have 

I>(5f )Va.(r) = XiKifi^a^ir), {p integer) (3) 

where X{Hii) is the eigenvalue of D{Si). Since electron 
density remains invariant under symmetry operations, 

|A(«i)P^/'aK(r)| = |VaK(r-)|, (4) 

we get X{hzi) = ex.p[ia{tvi)]. Now we impose periodic 
boundary conditions by making the group cyclic S^^ = 
1, and get a{hzi) = hzi = 27rmi/Mi, with integers Mi 
and mi G [0, Mi - 1]. 
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TABLE I: Selected examples on the approach usage. The coordinates in parentheses mean the sense of the symmetry operation. 



Operations 



Examples of usage & notes 



Si{z) rotation 

Si{z) rotation and S2{z) translation 

Si joined rotation(z) + translation(z) 

Si{x) and S2{y) two rotations around 
the same origin 

<Si(x), 82(1/) rotations around different 
origins 

Si translation(x)+reflection (yz) 

<Si(x), S2{y) translations, S^{xy) reflec- 
tion [+optional translation (x or y)] 

normal point group symmetries 



bend tubes, wires, ribbons [4j 
bend membranes, slabs 

twist nanotubes, wires, ribbons, DNA, simulate springs and coils p^tilG] 

simulate spherical symmetry; solid and liquid membranes, such as mono- and 
multilayer graphene and lipid bilayers. (<Si(x) and S2{y) commute approximately 
if rotation angles are small, and curvature can be treated as a local property.) 

simulating arbitrary Gaussian curvature, like saddle structures (approximate 
treatment, like above) 

structures with repeating units . . . ABBA . . ., using AB unit; many waves can be 
simulated using half the wavelength [17 

computational surface science; reflection doubles the surface slab thickness with 
half the number of atoms 

finite symmetric molecules and clusters 



By repeating the above steps for the remaining sym- 
metry operations, we obtain a revised Bloch's theorem: 
in a potential, that is invariant in symmetry operations 
tS'^, the energy eigenstates ifjan at r and at r' — S~'^r 
differ by a phase factor exp(zA^ • n), 

^(^")V^a.(r) = ^an^S-'^r) = exp(zK • n)^a.{r). (5) 

This implies that wave functions only in one unit cell — 
whatever its shape — determine the electronic structure 
of the whole, extended system. 

You may recognize that the theorem above is nothing 
but Bloch's theorem, merely written with unconventional 
symbols. Indeed, many things remain as usual. Energy 
eigenstates i/jan. automatically fulfill Eq.([5|, once written 
in a revised version of Bloch basis, 

^ n 

(6) 

where (Pjj,{r) are local orbitals and 1 = is the num- 
ber of unit cells. The Bloch basis gives Hamiltonian di- 
agonal in A^, 

(a^, /^\H\n\ v) = 5{n — n) ^ exp(— za^ • n)H^y(n) (7) 



with 



) = / ^;ir)H 



dV, 



(8) 



and similarly for overlap matrix elements. The total en- 
ergy expressions remain the same, we only use a set of 
K-points instead of /c-points (extra symmetries reduce 
the set). [18 Because forces are calculated as paramet- 
ric derivatives of the total energy, molecular dynamics 
works normally and energy is conserved; simulation cell 
dynamics, however, are different, and the concept of pres- 
sure needs redefinition. Finally, the theorem works with 
any electronic structure method, whether it is ab initio 
or not, whether it uses real-space grids or local orbitals 
(plane waves are tricky), or whether the approach is nu- 
merical or analytical. 



Some things, however, do change in the revised Bloch's 
theorem. For bulk the periodic boundary condition is 
an approximation, whereas here some symmetries may 
form cyclic groups in reality, as in benzene. For cyclic 
groups the K-point sampling is more restricted; in the 
above example of an achiral carbon nanotube the trans- 
lational component i^i can be freely sampled between 
[— 7r,7r] (because for Si periodicity is an approximation), 
but the rotational component accepts only the discrete 
values hi2 = 27rm/M, m = 0, 1, . . . , M — 1. Group mul- 
tiplication tables for Si that have identities like S\^ = Sj' 
make the sampling of the components of k coupled, for 
then we must have = IjKj + 27rm (/i,m integers). 
The connection between k- and /c-points is as follows. If 
Si is a translation along L^, then Siipan should give the 
same phase as T{L^)il^ak^ if i^an and il^ak are the same 
physical states. Hence exp(— mi) = exp(— z/c ■ L^)^ and, 
in general, tii = Y^j^i ^ ~ 1, 2, 3. 

The formalism gives surprises, too. An atom can per- 
form work on itself. This is because the total force on 
an atom exerted by its own periodic images — if rotations 
are involved — may differ from zero. Furthermore, a force 
Fji on atom / exerted by atom J may not be the coun- 
terforce to the force on atom J exerted by atom /, that is 
F ji 7^ —Fij] Newton's third law appears invalid. These 
unorthodoxies are not bugs; remember that we simulate 
the whole extended system, and an atom / in the prim- 
itive unit cell is different from the atom / in a different 
unit cell — an artifact of atom indexing. Finally, note that 
if Si contains rotations, also local orbitals rotate; this is 
implicit in the operation D{Si)(p^{r) in Eq.|6|. 

We implemented this formalism using localDasis in the 
density- functional tight-binding software hotbit p!9[ [20]. 
and tested it with many finite and extended structures. 
We omit the details of the implementation here, and just 
comment on three things. First, the standard methods 
of electrostatics, like Ewald summation, are invalid since 
flexibility is required; we chose to use multipoles as they 
easily lend themselves for rotations and reflections. Sec- 
ond, the implementation can be done so that only the 
mappings 
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are needed to build new symmetries; this requires just 
a couple of lines new code. Implementation generally 
is not hard, but it may be nontrivial for codes already 
build upon translational symmetry. Third, implementa- 
tion has a negligible computational overhead compared 
to translational symmetry (see Table Certain ma- 
nipulations take more time, but the most CPU-intensive 
parts remain as usual. 

So far our discussion has been abstract, but what can 
we do with the formalism in practice? While it may 
seem that we require a lot of symmetries, the main point 
of this Letter quite the opposite: we require less symme- 
tries than before. Formalism enables simulating distorted 
materials, but also reduces computational costs for cer- 
tain simulations. Selected examples of usage are shown in 
Table |Tj For example, one cost-reduction area is surface 
science, where less atoms are needed to simulate thick 
surface slabs. We believe more application areas can be 
discovered, once the new concepts are mastered. 

Now we leave the general discussion, and give one prac- 
tical example of usage: we investigate twisted armchair 
graphene nanoribbons (AGNRs).[21 We choose this ex- 
ample for the existing literature, but also for the pos- 
sibility to illustrate operations beyond standard chiral 
symmetry. 

Figure [l^ shows a piece of an infinitely long 20-AGNR, 
with a twist x = 360 deg/290.7 A= 1.2 deg/A, within 
one conventional unit cell of 2720 atoms. The minimal 
unit cell, enabled by the new formalism, in turn, has 
20 atoms (atoms A in Fig. [T]3), and is accompanied by 
two symmetry operations: Si is L = 4.2 A translation, 
followed by L • X = 5.04 degree rotation, and S2 is L/2 = 
2.1 A translation, followed by 182.52 degree rotation (= 
180 deg+x ' The whole system can be built from 

one unit ceh by S]^'Sp with mi = 0,±1,±2, ... and 
m2 = 0, 1. Note that S2 = Si and hence = 1^1 / 2 -\- 
7rm2, while hci is freely sampled. (While this was our 
choice for the symmetry operations, another, and equally 
sufficient, choice would have been to use only ^2 with 
m2 =0,±1,±2,....)[22] 

Table [TTl shows the wall-clock times for selected simu- 
lations. The simulations with x = show that the new 
formalism has no computational overhead compared to 
translation. The simulations with x = 1-2 deg/A show 
that finite twist affects simulation times with neither 
minimal nor chiral cells. The translational cell was too 
large for direct simulation, and the timing was estimated 
from the scaling law time~ (system size)^. Note that, by 
decreasing x? the translational cell size — along with its 
timing — could easily be grown indefinitely. Point here is 
that twists even smaller than 1.2 deg/A will be required 
to investigate the relevant physics of a 20-AGNR. Con- 
ventional quantum-mechanical simulation is practically 
impossible. We also remark that, given compatible Ap- 
point samplings, energies and forces from different types 
of cells are the same within ffoating-point precision. 

Figure ^ shows AGNRs' energies as a function of 
twist. For ribbons wider than 12 A the energy is at 
minimum with non-zero X = Xo; ribbons twist sponta- 
neously. This confirms earlier predictions by classical 
potentials (using thousands of atoms) [5l |9] and finite el- 
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FIG. 1: (Color online) (a) Conventional unit cell of 20-AGNR 
with a twist x = 1-2 deg/A (2720 atoms), (b) Minimal unit 
cell (atoms A), with illustrations of Si and S2- Si is the chiral 
operation and S2 transforms A^B. (c) Energies for AGNRs 
as a function of x for selected ribbon widths. Minimum of 
the curve gives the spontaneous twist xo- We used a 20 x 2 
K;-point mesh. Inset: xo as a function of AGNR width. 



TABLE II: Timings for 20-AGNR with different twists. Time 
is the wall-clock time per molecular dynamics step, calculated 
with a standard present-day desktop computer. Timing 
obtained from scaling law. 



unit cell 



atoms X (deg/A) time (s) 



minimal (like A in Fig. lb 
chiral (like A+B in Fig.^ 
translational (like chiralj 
minimal (A in Fig. Ip) 
chiral (A+B in Fig. it) 
translational (Fig, [l 1) 



20 





1.51 


40 





2.75 


40 





2.75 


20 


1.2 


1.51 


40 


1.2 


2.75 


2720 


1.2 


9.26 X 10 



ement modeling [3 . The physical reason for twisting is 
the compressive edge stress that elongates edges with re- 
spect to ribbon's center. [5l [H [23l |Mj The stress we get 
(1.7 eV/A) agrees well with the stress (1.5 eV/A) from 
previous density-functional calculations. [9^ For wide rib- 
bons we get scaling xo ^ 22 deg/width, and the differ- 
ence to classical scaling xo 18 deg/width of Ref. [5] 
comes mainly from quantum mechanics: the edge stress 
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FIG. 2: (Color online) (a) Energy gap of a 10.8 A wide 10- 
AGNR as a function of x- The 7r-only gap is obtained by 
treating the unhybridized 2p;z -electrons as s-electrons. Inset: 
energy gaps as a function of AGNR width for x = 0. (b) Band 
structures for flat and twisted ribbons, plotted as a function 
of Ki. The symmetry operation S2 brings another dimension, 
to band structure plots. Calculation does not include spin. 



resides not only at the edge, but extends more into rib- 
bon's center — a feature hard to reproduce by classical 
potentials. [9 This is also why we have no spontaneous 
twist for narrow ribbons. Ribbons ^ 11 A wide have 
nearly zero torsion constant, and could be used in ultra- 
sensitive torsion balances. Ultimately, very wide ribbons 
should show bifurcation into flat ribbons with ripples at 
the edges, but we won't discuss that here.^l£] 

As argued in the abstract, the electronic properties 
ought to be investigated together with mechanical prop- 
erties; this requires quantum mechanics. It is known that 
AGNRs have a gap due to the confinement of the finite 
width, and our gaps (inset in Fig. |2^) agree well with 
density- functional calculations of Reir21, But what hap- 
pens to electronic structure when ribbons get twisted? 
Fig. [2^ shows that twisting changes 10-AGNR's gap very 
little — this is generic for all AGNRs. The gap from tt- 
electrons alone shows further that sp-rehybridization is 
negligible. Even the band structures of flat and twisted 
ribbons (Fig. [2b) are nearly identical. This suggests that, 
contrary to CNTs[25 , the electronic properties of GNRs 



are remarkably robust against twisting. 

To conclude, we hope to have illustrated how mod- 
est revision of Bloch's theorem enables versatile material 
distortions with quantum mechanics included, both nu- 
merically and analytically. However, excess emphasis on 
quantum mechanics causes undue discrimination of clas- 
sical methods — the formalism works equally with classi- 
cal force fields, finite element methods, or coarse-grained 
simulations, and equally when applied to, say, liquid- 
phase cell membranes, fluid flow through bent pipes, or 
electron transport. 
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